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The Birkhoff Ergodic Theorem concludes that time averages, that is, Birkhoff averages, T,^ =1 f(x n )/N 
of a function / along an ergodic trajectory ( x n ) of a function T converges to the space average f fd/i , 
where fi is the unique invariant probability measure. Convergence of the time average to the space 
average is slow. We introduce a modified average of f(x n ) by giving very small weights to the “end” 
terms when n is near 0 or N. When (x n ) is a trajectory on a quasiperiodic toms and / and T are C°°, 
we show that our weighted Birkhoff averages converge “super” fast to / /d/i, i.e. with error smaller 
than every polynomial of 1/N. Our goal is to show that our weighted Birkhoff average is a powerful 
computational tool, and this paper illustrates its use for several examples where the quasiperiodic set 
is one or two dimensional. In particular, we compute rotation numbers and conjugacies (i.e. changes 
of variables) and their Fourier series, often with 30-digit precision. 


Introduction 


Quasiperiodicity is a key type of observed dynamical behavior in a diverse set of applications. Tori with 
(^jasiperiodic motion persist for small perturbations by the Kolmogorov-Arnold-Moser theory, but such 
<J>e|havior is also observed for non-conservative systems well beyond this restricted regime. We believe 
gat quasiperiodicity is one of only three types of dynamical behaviors occurring in basic sets of typical 
©stems. See [lj for the statement of our formal conjecture of this basic set triumvirate. For example, 
^asiperiodicity occurs in a system of weakly coupled oscillators, in which there is an invariant smooth 
tracting torus in phase space with behavior described exclusively by the phase angles of rotation of 
system. Indeed, it is the property of the motion being described using only a set of phase angles 
always characterizes quasiperiodic behavior. In a now classical set of papers, Newhouse, Ruelle, and 
Tokens demonstrated a route to chaos through a region with quasiperiodic behavior, causing a surge in 
study of the motion PJ. There is active current interest in development of a systematic numerical and 
tfhjeoretical approach to bifurcation theory for quasiperiodic systems. 

Our goal in this paper is to present a fast numerical method for the fast calculation of the limit of 
Birkhoff averages in quasiperiodic systems, allowing us to compute various key quantities. If / is integrable 
and the dynamical system is ergodic on the set in which the trajectory lives, then the Birkhoff Ergodic 
Theorem asserts that the Birkhoff average T,% = 1 f(x n )/N of a function / along an ergodic trajectory 
(x n ) converges to the space average f fd/j, as N -> oo for //-almost every xo, where // is the unique 
invariant probability measure. We develop a numerical method for calculating the limit of such averages, 
where instead of weighting the terms f(x n ) in the average equally, we weight the early and late terms 
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of the set {1,... ,N} much less than the terms with n ~ N/2 in the middle. That is, rather than using 
the equal weighting (1 /N) in the Birkhoff average, we use a weighting function w(n/N ), which will 
primarily be the following well known C°° function that we will call the exponential weighting function, 
w exp (t) = exp(l/(i(t- 1)). In a companion paper [3], it is rigorously shown that for quasiperiodic systems 
and C°° /, this weighting function leads to super convergence with respect to N , meaning faster than 
any polynomial in N~ l . This super convergence arises from the fact that we are taking advantage of 
the quasiperiodic nature of the map or flow. In particular, our method uses the underlying structure 
of a quasiperiodic system, and would not give improved convergence results for chaotic systems. We 
demonstrate the method and its convergence rate by computing rotation numbers, conjugacies, and their 
Fourier series in dimensions one and two. We will refer to a ID quasiperiodic curve as a (topological) 
circle. 

Other authors have considered related numerical methods before, in particular 00, which we will 
compare to when we introduce our method. See also 0 m El U EDI HHH21H3UH US] . 

Our paper proceeds as follows: In Section [2j we give the formal definition of quasiperiodicity, rotation, 
and the conjugacy map. In Section |3.1[ we describe our numerical method in detail. We illustrate our 
method for a series of four examples, including an example of a two-dimensionally quasiperiodic map. In 
all cases, we get fast convergence and are in most cases able to give results with thirty digits of precision. 
For convenience of the reader, have summarized our numerical findings in Table [lj 

We start by describing our results for a key example of quasiperiodicity: the (circular) restricted three- 
body problem. This is an idealized model of the motion of the planet, a large moon, and a spacecraft 
governed by Newtonian mechanics, in a model studied by Poincare mm In particular, we consider 
a planar three-body problem consisting of two massive bodies (“planet” and “moon”) moving in circles 
about their center of mass and a third body (“spacecraft”) whose mass is infinitesimal, having no effect 
on the dynamics of the other two. 

We assume that the moon has mass p and the planet mass is 1 - // where p = 0.1. and writing equations 
in rotating coordinates around the center of mass. Thus the planet remains fixed at (-0.1,0), and the moon 
is fixed at (0.9,0). In these coordinates, the satellite’s location and velocity are given by the generalized 
position vector ( (p , q 2 ) and generalized velocity vector (pi,p 2 )- The equations of motion are as follows (see 

inn- 


dqi/dt = Pi + q 2 , 
dq 2 /dt = p 2 ~qi, 

dpi/dt = p 2 — p{qi - 1 + p)d^ oon - (1 - p) (q 1 + p)d~f anet , 


( 1 ) 


dp 2 /dt = -pi - pq 2 d. 


-3 

moon 


(1 - p)q 2 d, 3 


'planet ’ 


where 


dmoon ~ (((Zi 1 a0 + 0.2 ) ctnd d p i anet — ((<Ai + a0 + ^ 2 ) 


2\0.5 


The following function H is a Hamiltonian for this system 


H = [{pi +pl)/2] + [q 2 pi - qip 2 ] 


[P dpianet + C 


I- 1 ) dj 00 J. 


( 2 ) 


The terms in the square brackets are resp. the kinetic energy, angular moment and the angular potential. 
For fixed H, Poincare reduced this problem to the study of the Poincare return map for a fixed value 
of H, only considering a discrete trajectory of the values of (qi,pi) on the section q 2 = 0 and ^ > 0. 
Thus we consider a map in two dimensions rather than a flow in four dimensions. Figure [2] shows one 
possible motion of the spacecraft for the full flow. The orbit is spiraling on a torus. The black circle shows 
the corresponding trajectory on the Poincare return map. Fig. [l] shows the Poincare return map for the 
spacecraft for a variety of starting points. A variety of orbits are shown, most of which are quasiperiodic 
invariant circles. An exception is A-trajectory in Fig. [IJa), which is an invariant recurrent set consisting 
of 42 circles. Each circle is an invariant quasiperiodic circle under the 42-nd iterate of the map. Using our 
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Example 

Equation 

Rotation number (s) 

Relatec 


7 

igures 

Restricted three-body problem 

1 


0.063961728757453097164077919081024 

T] 

3 


2| 

Standard map 


2 

0.12055272197375513300298164369839 

T 


5 


Forced van der Pol oscillator, F = 5 

1 

3 

0.29206126329199589285577578718959 


H 

6 



Forced van der Pol oscillator, F = 15 

1 

3 

0.37553441113144010884908928083318 
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Forced van der Pol oscillator, F = 25 

1 

3 

0.56235370092685056634419221336154 
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Two-dimensional torus 

1 

4 

Pi , p 2 in Table [ 2 ] 

_ E 

0 

9 

1 _ 


Table 1: Summary of our numerical calculations. 


numerical method for Birkhoff averages, in the itemized list given below we summarize our results for the 
quasiperiodic orbit labeled B\. 

The error in the quasiperiodicity computations using weighted Birkhoff averages decreases exponen¬ 
tially in the number of iterates N (see Fig. [3^) . This speed of convergence means the accuracy of our 
solutions is the limit of numeric precision. In particular, we have computed trajectories for the Poincare 
return map using an 8-th order Runge-Kutta method with time step 10 -5 , in quadruple precision, mean¬ 
ing our results given below are computed up to thirty digits of accuracy. Section [2] formally defines the 
computed values given in this list. 

1. The rotation number is given in Table [1} computed to 30 digits of accuracy. 

2. We can compute the Fourier series of up to 200 terms. There is a conjugacy map h between the first 
return map and a rigid rotation on the circle. Evaluating the Fourier series allows us to reconstruct 
the conjugacy map (c/. Fig. [djr). 

3. The exponential decay of the coefficients in the Fourier series is a strong indication of the analyticity 
of conjugacy function (c/. Fig. [3}> ). 
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( c ) 

Figure 1: Poincare-return map for the restricted three-body problem. All three parts show a 
projection to the q\ -pi plane. Figures (a) and (b) show various quasiperiodic trajectories on the Poincare 
section q 2 = 0, of the restricted three-body problem ([!]). Note that the planet is fixed at the point (-0.1,0) 
and the planet at (0.9,0). Thus some trajectories orbit both the planet-moon system and some only orbit 
the planet or the moon. Each time the flow hits q 2 = 0 and ^ > 0, we plot (qi,pi). Each trajectory 
shown is a (topological) circle with quasiperiodic motion. The energy H for all the circles shown in the 
figures, including B\ is the same and H « -2.63. Part (c) shows in white all the initial points (qi,pi) on 
the Poincare surface for which there exists a p 2 so that the Hamiltonian H at {(p , q 2 = 0,pi,p2) is the same 
as the one in parts (a) and (b). Part (c) also shows the trajectory which corresponds to the circle B\ in 
(a) and (b). 
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Figure 2: Torus flow for the restricted three-body problem. All four views are of the same two- 
dimensional quasiperiodic torus lying in R 4 . Each picture consists of the same trajectory spiralling densely 
on this torus. This trajectory is the solution of ([!]), shown as curve B\ in Fig. [lj We require four different 
views of this torus because the embedding into three dimensions gives a highly non-intuitive images. The 
black circle is the set of values of the Poincare return map with q 2 = 0 for this flow torus. 
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Figure 3: Quasiperiodicity for the restricted three-body problem. For the quasiperiodic circle B\ 
in Fig. [lj Part(a) shows how the invariant circle is parameterized by coordinates 0 € S l = [0,1). Part (b) 
depicts the periodic part g{6) of the conjugacy between the quasiperiodic behavior and rigid rotation by 
p. See 0 for a description of g(0). Part (c) shows the norm of the Fourier coefficients of the conjugacy 
as a function of index. This exponential decay indicates that the conjugacy function is analytic, up to 
numerical precision. Part (d) shows the convergence rate of the error in the rotation number p m as a 
function of the number of iterates N. The “error” is the difference | p^ - Pn* |, where N* is large enough 
so that pm appears to have converged. The step size used for the Runge-Kutta (RK8) scheme is 10 -5 . 
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2 Quasiperiodicity 

In the introduction, we described quasiperiodic motion as motion that could be fully understood through 
a set of angles of rotation. We now formalize that idea in the following definition for quasiperiodicity. 

Quasiperiodicity For a dimension d > 1, let p = (pi,p 2 , • • • ,Pd) he a rotation vector such that all pk 
are irrational and irrationally related (defined below in (|8j) ). Then following map is known as a rigid 
irrational rotation: 


Tp(0) = 6 + p (mod 1), where 0 < 9k < 1, k = 1,2,... ,n. (3) 

A rigid rotation the simplest, albeit least interesting example of a map with quasiperiodicity. Since Tp 
gives the same values on opposite sides of the unit cube, we identify the sides and refer to the domain of 
the rigid rotation as a circle in one dimension and an n-torus in n > 1 dimension. From now on we will 
refer to the circle as a 1-torus. We define a general map T to be quasiperiodic if either T or some iterate 
T k is topologically conjugate to a rigid rotation. (We will assume k = 1 in the rest of this description.) 
That is, a map T is quasiperiodic if there is a rigid rotation map Tp and an invertible conjugacy map h 
such that 

T(h(0)) = h(T ls (0)). 

A flow has quasiperiodic behavior if its associated first return Poincare map has quasiperiodic behavior. 

Since we are focused on maps, we restate this definition in terms of iterates. Let ( x n ) be the forward 
orbit under T, and ( 9 n ) the forward orbit under Tp. That is, x n+1 = T(x n ) and 9 n+ \ = 9 n + p (mod 1). 
Then as long as 0o = h($o) = h( 0), 


(j) n = h(9 n ) for all n = 0,1, 2,3,.... (4) 

The map h maps a torus to the domain of T, meaning that the domain of T is a topological circle or torus. 
Since h is a homeomorphism, the map 

g(f>) :=h(0)-f) (5) 

is periodic. 

Thus a map T restricted to an invariant topological circle C is quasiperiodic if there is an continuous 
invertible map h from C to the true circle which maps T to a rigid irrational rotation map. 

For an invertible map F to be quasiperiodic on a topological circle C. it is necessary and sufficient 
that F has a dense trajectory, as shown in ng. In general, a circle homeomorphism without periodic 
points may not be quasiperiodic. However, if we assume that the map F and the topological circle C 
are twice continuously differentiable, then Denjoy [19] showed that these conditions are both necessary 
and sufficient. Furthermore, clearly any rigid irrational rotation map is a real analytic map, but even 
if we assume that a quasiperiodic function is analytic, Arnold showed that the conjugacy map h may 
only be continuous for some atypical rotation number. However, Herman (see [20]) proved that for circle 
homeomorphisms, most conjugacy maps h are analytic. 
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3 Our WBjv+ method and its applications 


3.1 Rotation number 


Here is our method for calculating the limit of the Birkhoff average lim^oo S n=if( x n)/N = f f dp along 
an ergodic trajectory (or first return map) ( x n ). Let the weighting function w be defined as 


Wexp(t) ■= w(t) := 



for t € (0,1) 
for t i (0,1). 


( 6 ) 


Let T be quasiperiodic on some set X 0 , with X l} a topological torus. Henceforth, X 0 = T d will denote a 
d-dimensional topological torus. For d = 1, T d is a topological circle. For a continuous function / and a 
C°° quasiperiodic map T on T d , let x n € T d be such that x n = T(x n - 1 ) for all n > 1. We define a Weighted 
Birkhoff (WBjv) average of / as 


N -1 

WBa r(f)(x) := Y, Wn,Nf(xn), where ui n)N 

n =0 

In a companion paper [3] , it is shown that if / and T are C °°, 
is a constant C m > 0 such that 


w(n/N ) 

£?.XiW 


then for every positive integer m there 


WB N (f)(x) - J fdn 


< C m N m for every m. 


We refer to the above as super (polynomial) convergence or exponential convergence. 

Our numerical method converges fast enough to allow us to obtain high precision values for f fdfi 
with relatively low computational cost. In particular, by appropriate choices for /, we obtain rotation 
numbers and sometimes the Lyapunov exponents, Fourier coefficients, and conjugacy reconstructions for 
quasiperiodic maps and flows potentially in any finite dimension, though here we work mainly in dimension 
one, with one example in dimension two. We now show how to apply this general method to computation 
of specific quantities. We observe that N must be larger for T 2 than for T 1 , to get a 30-digit accuracy. 

Calculation of rotation vectors. Recall that every quasiperiodic map is conjugate to a pure 
rotation T p by a vector p = (pi...., p n ). In our discussion, we assume that the components of p are 
irrational and irrationally related. When we say that pi,...,pd are irrationally related, we mean that 
there are no integers ki ,..., kd and Ay for which, 


A’iPi + k2P2 + —t- kdPd + ko 


( 8 ) 


unless all ki are zero. Since we are only provided with the map T and not with the map R, we need a way 
of recovering this rotation vector using only a forward trajectory of the map T. The following formula 
gives a way of calculating the rotation vector p directly from the original map. The rotation vector given 
in the formula is valid well beyond the case of a quasiperiodic function. For example, the rotation vector 
can certainly be rational in one or all of its components. 

Rotation vectors. Let q(x) = x (mod 1) and define q : ]R d -> T d to be the projection map modulo Z in 
each coordinate. For example, in one dimension the circle [0,1] with ends identified is the image of the 
real M under the map q. For any homeomorphism T : T d ->■ T d . let T : M. d ->■ be the lift of T. That is, 

q o T = (f) o q. 

Since T is a homeomorphism, the map T(x) -x must be periodic in x with period 1. The rotation vector 
of a homeomorphism T : T d -»• T d starting at lift point z is given by 








(a) 


(b) 


Figure 4: The standard map. A variety of orbits from different initial conditions in the standard map 
Si defined in (12) are plotted on the left. We can see both chaos (in pink) and quasiperiodic orbits under 
this map. A single topological circle with quasiperiodic behavior is plotted on the right. The orbit has 
initial conditions (x,y) » (-0.607,2.01).That is, if we restrict the map to this invariant circle, then it is 
topologically conjugate to a rigid irrational rotation. 


Note that possible rotation numbers (p i. /q) of' a quasiperiodic map on a torus T d depend on the coordinates 
chosen for T and z. For example, for domain T 2 , the set of pairs (pi, P 2 ) possible for a fixed T is dense in 
the two-torus. See [3j. 

Based on the above definition, the standard method of computing p from an orbit (y^) of length n is 
to calculate the average 

P* 4 Z(^(w)-2/n)- (10) 

iV n-1 


This method converges slowly at best, with order of only 1/N |4j. However, since Eq. 10 can be written 
as a Birkhoff average by writing f(y n ) = T(y n ) - we can apply our method to this function. That is, 
let (yk)k=o be a forward orbit for T. Our approximation of p is given by the weighted average of /, 


N -1 

WB N (y n+ l - IJn) ■= y ™n,N(y n + 1 - y n ) P- (H) 

n =0 

Rotation number for the restricted three-body problem. In the introduction, we stated our 
results, including the rotation number for a quasiperiodic orbit for the restricted three-body problem. 
Fig. [3](b) shows the convergence rate of the calculation for this rotation number. . 

Rotation number for the standard map. The standard map is an area preserving map on the 
two-dimensional torus, often studied as a typical example of analytic twist maps (see EH). It is defined 
as follow^] 

/*W x + v \ {mod2n) (12) 

\ y j \ y + sm xj 

The left-hand panel of Fig. |4] shows the trajectories starting at a variety of different initial conditions 
plotted in different colors. The shaded set is a large invariant chaotic set with chaotic behavior, but 


*The standard map generally depends on a parameter a , but we only consider the case a = 1. 


9 














Figure 5: The Standard map conjugacy. Fig. (a) shows the shows the change of variables which 
converts the motion on the circle in Fig. & into a pure rotation. Fig. (b) shows the decay of the Fourier 
coefficients. Since the conjugacy is an odd function, the odd-numbered Fourier sine and cosine terms are 
zero and therefore have been omitted from the picture. The decay of the Fourier terms can be bounded 
from above be an exponential decay, which suggests that the conjugacy is analytic. 


many other invariant sets consist of one or more topological circles, on which the system has quasiperiodic 
behavior. For example, initial conditions (ir, 1.65) leads to chaos while (n, 1.5) leads to a quasiperiodic 
trajectory. As is clearly the case here, one-dimensional quasiperiodic sets often occur in families for non¬ 
linear processes, structured like the rings of an onion. There are typically narrow bands of chaos between 
quasiperiodic onion rings. Usually these inner rings are differentiable images of the n-torus, Yamaguchi 
and Tanikawa EH and Chow et. al. in [22] showed that the outermost limit (the onion’s skin, to continue 
the analogy) will still be a torus, but may not be differentiable. We have computed the rotation number 
for the standard map orbit shown in the right panel Fig. [4] using quadruple precision. Its 30 digit value is 
given in Table [T} 

The forced Van der Pol oscillator. Fig. [6] shows orbits for the time-27r/0.83 map of the following 
periodically forced Van der Pol oscillator with nonlinear damping [23] 

- 0.2 (l - x 2 ) + 20x 3 = Fsin (0.83 1 ), (13) 

H/L 

with a variety of choices of F. While the innermost orbit shown is a chaotic attractor, the outer orbits are 
topological circles with quasiperiodic behavioijt} Our computed the rotation number for the three orbits 
F = 15,25,35 are given in Table [l] 

A two-dimensional torus map. We now introduce an example of a two-dimensional quasiperiodic 
torus map on T 2 . This is a two-dimensional version of Arnold’s family of ID maps (see [24]), originally 
introduced in the following two papers [251126]. The map is given by (Ti,T 2 ) where 


Ti(x,y) 

T-z(x,y) 


X + ^1 + ^Pi(x,2/)) 
Z7T 

y + uj 2 + ^-P 2 (x,y) 


(mod 1) 
(mod 1), 


As with the standard map, we have specified all non-essential parameters rather than stating the most general form of 
the Van der Pol equation. 
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Figure 6: Forced van der Pol oscillator. Part (a): Attracting orbits for a number of different forcing 
values F for the stroboscopic map of the van der Pol flow given in (13). The plot depicts points (A, Y) = 

(x(tk),x'(tk)), where tj~ = 2£)7r/0.83,fc = 0,1)2,_The chaotic orbit lying inside the cycles corresponds 

to F = 45.0. There are stable quasiperiodic orbits shown as circles, which from outermost to innermost 
correspond to F = 5.0, 15.0, 25.0 and 35.0 respectively. Part (b): the periodic part g(6) of the conjugacy 
(Eq. i to a pure rotation, for F = 5, 15 and 25. Part (c): the second figure gives a demonstration of 
the analyticity of the conjugacy to a pure rotation by studying the decay of the norm of the k -th Fourier 
coefficients with k, up to the resolution of the numerics. 


and Pi(x,y),i = 1,2 are periodic functions with period one in both variables, defined by: 

4 

Pi(x, y) = Yj a i,j sin( 27 TQ! i j), with a itj = rjx + sjy + b itj . 

3 = 1 

The values of all coefficients are given in Table [2j This choice of this function is based on [251 [26]. Both 
papers use the same form of equation, though the constants are close to but not precisely the same as 
the ones used previously. This is fitting with the point of view advocated in by these papers, in that the 
constants should be randomly chosen. Since we are using higher precision, we have chosen constants that 
are irrational to the level of our precision. The forward orbit is dense on the torus, and the map is a 
nonlinear which exhibits two-dimensional quasiperiodic behavior. 

Fig. [T^l depicts iterates of the orbit, indicating that it is dense in the torus. In order to verify that this 
map is indeed quasiperiodic, we used a Birkhoff average to compute the two Lyapunov exponents to verify 
that both are zero (c/. Fig. &) In terms of method, this is just a matter of changing the function / used 
in WB ; \? in Eq. [7] Likewise, finding rotation numbers in two dimensions uses the same technique as in 
the one-dimensional case (c/. Fig. [?[■) . In all of our calculations, the computation is significantly longer 
than in one dimension in order to get the same accuracy, perhaps because in two dimensions, coverage of 
dense orbit varies like the square of the side length of the domain. 

Convergence rate. In order to explain why the convergence of our method is so good, we introduce 
four different possible values for the weighting function w, depicted in Fig.jSjr, and compare the convergence 
results for computing the rotation number for each of these choices of w. 


^equal (0 

= 1 (Birkhoff’s choice) 

(14) 

^ quadi )0 

= 

(15) 

^(sin 2 )(0 

= sin 2 (7rf) 

(16) 

^exp(^ ) 

= exp (-l/(t(l -£)))• 

(17) 
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Coefficient 

Value 

€ 

0.4234823 

kfi 

0.71151134457776362264681206697006238 

0 J 2 

0.87735009811261456100917086672849971 

a i ,3 

(-0.268,-0.9106,0.3,-0.04) 

a 2 ,j 

(0.08,-0.56,0.947,-0.4003) 

1—1 

(0.985,0504,0.947,0.2334) 


(0.99,0.33,0.29,0.155) 

r j 

(1,0,1,0) 

S j 

(0,1,1,-1) 

Computed pi 

0.718053759982066107095244936117 

Computed P 2 

0.885304666596099792113366824157 


Table 2: Coefficients for the torus map. All values are used in quadruple precision, but in this table 
the repeated zeros on the end of the number are suppressed. 


If we compute with the first choice of w, we recover the a truncated series in the definition of the Birkhoff 
average. To estimate the error, we expect the difference /(a^jv+i) _ /(^v) to be of order one, implying that 

WBw equal ,N +1 - WB w equal ,N ~ 1/N. 

The choice of a particular starting point also creates a similar uncertainty of order 1/N. For all but the 
first choice, w is always positive between 0 and 1 but vanishes as t approaches 0 and 1. In addition, going 
down the list, increasing number of derivatives of w vanish for t -> 0 and t -> 1, with all derivatives of w exp 
vanishing at 0 and 1. We thus expect the effect of the starting and endpoints to decay at the same rate as 
this number of vanishing derivatives. Indeed, we find that u) qU ad corresponds approximately to order 1/N 2 
convergence, u>( sin 2 ) to 1/N 3 convergence, and w exp to convergence faster than any polynomial. Figs. [8|) 
and [7bc show this effect. 

Related methods. See HI El for references to earlier methods for computing rotation numbers. 
In [H 5], A. Luque and J. Villanueva develop fast methods for obtaining rotation numbers for analytic 
functions on a quasiperiodic torus, sometimes with quasiperiodic forcing with several rotation numbers. 
The paper [5] examines a smooth function f on a quasiperiodic torus. Let f n denote the value of f at 
the n-th trajectory point. From this sequence they can obtain the rotation number with error satisfying 
|error| < C p N~ p for any p where C p is a constant. The method of computation depends on p and as p 
increases the computational complexity increases for fixed N. If T(p,N ) is their computation time, it 
appears that T(p,N)/N -»• oo as p -»• oo. In comparison, computation time for our weighted Birkhoff 
average is simply proportional to N since it requires a sum of N numbers. The paper gives one figure (Fig 
6) from which the rate of convergence can be computed, namely a restricted three-body problem. Their 
rotation-rate error is proportional to A -3 - 5 and is » 10 -18 at N = 2 21 . 

Several variants of the Newton’s method have been employed to determine quasiperiodic trajectories 
in different settings. In m the monodromy variant of Newtons method was applied to locate periodic 
or quasi-periodic relative satellite motion. A PDE-based approach was taken in [28], where the authors 
defined an invariance equation which involves partial derivatives. The invariant tori are then computed 
using finite element methods. See also Chapter 2, [28] for more references on the numerical computation 
of invariant tori. 
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3.2 Fourier coefficients and conjugacy reconstruction 

For a quasiperiodic curve as shown in Fig. [3}i, there are two approaches to representing the curve. Firstly, 
we can write the coordinates (X,Y) as a function of 9 e S 1 , or secondly, we can reduce the dimension and 
represent the points on the curve by an angle <f> e S 1 , that is, cf)(X(9),Y ( 9 ), which is also h(9) = 9 + g{9). 
We have shown g in Fig [3]b and the exponential decay of the norm of the Fourier coefficients in Fig. [ 3 J 3 . 
To limit the number of graphs in this paper, we have only created the Fourier series for the periodic part 
9(0) ' ‘ . 

Given a continuous periodic map / : S 1 -*■ M, where S' 1 is a circle (or one-torus), the Fourier sine and 
cosine representation of / is the following. 


For every t e S 1 , f(t) = l 

Z k= 1 


(18) 


k =0 


where the coefficients b k and c k are given by the following formulas. 


bk = 2 f f (9) cos(2kn9) d9, 

JOiS 1 

c k = 2 J f(9) sin(2fc7r$) d9. 


(19) 


( 20 ) 


To be able to use the fast Fourier transform, 2 M equally spaced points on the circle are required. If 
we only have access to an ergodic orbit (x n ) on the circle, then we cannot use the fast Fourier transform 
as we only have the function values f(x n ) along a quasiperiodic trajectory, and a rotation number p. So 
instead, we obtain these coefficients using a weighted Birkhoff average on a trajectory ( x n ) by applying 
the functional WB^jv- For k = 0, we find ao by applying WB to the function 1. For k > 0, we find b k and 
c k as follows. 


N 


b k = WB WtN (f(9)cos(2kn9)) = Y f( x n) cos(2kimp)w ntN . 

n -0 


N 


c k = WB„,»(/(6) Sin(2*:7r#)) = y f(x n ) sin(2kmp)w UiN . 


( 21 ) 


( 22 ) 


n -0 


By specifying that 6q = 0, our computation of rotation number p provides all iterates : 9 n = np. Using the 
Fourier coefficients, we can thus reconstruct the periodic part of the change of variables function g (see 
Eq. [5]). This is depicted for the restricted three-body problem in Fig. [3j for the standard map in Fig. [5| 
for the forced van der Pol equation in Fig. 6j In all three one-dimensional cases, we depict \Jb 2 k + c 2 k as a 
function of k. Our main observation is that the Fourier coefficients decay exponentially; that is, for some 
positive numbers a and f3, in dimension one, the Fourier coefficients b k and c k satisfy 


V\h | 2 + |c fc | 2 < ae for all k € Z. 


(23) 


This is characteristic of analytic functions. We therefore state that all of the conjugacy functions that 
we computed in our examples are effectively analytic, “effectively” meaning within the precision of our 
quadruple precision numerics. 

In two dimensions, the computation of Fourier coefficients is similar, but instead of only having one 
set of cosine and sine functions, for each (j, k), we have two linearly independent sets of complex-valued 
functions, where i = \f-l\ 

e i(jx+ky) anc j e i(jx-ky)' 

We define a hk and b hk to be the complex-valued coefficients corresponding to each of these functions. The 
reconstructed conjugacy function and decay of coefficients for the two-dimensional torus is depicted in 
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Fig. [ 9 J The decay of coefficients shows \fj 2 + k 2 on the horizontal axis, and \b :) j~\ and \cj )k \ on the vertical 
axis, where both of these coefficients are complex, meaning that |-| represents the modulus. Again here, the 
coefficients decay exponentially, though the decay of coefficients is considerable slower in two dimensions 
due to the added dimension. The data looks quit e a lot more crowded in this case, since there are many 
different values of (j, k) such that the values of \Jj 2 + k 2 are identical or very close. In addition, the two 
sets of coefficients bj :k and c ]k generally converge at different exponential rates. This is why there is a 
strange looking set of two different clouds of data in Fig. [9]b. While more information on the difference 
between these coefficients is gained by interactively viewing the data in three dimensions, we have not 
been able to find a satisfactory static flat projection of this data. We feel that in a still image, the data 
cloud shown conveys the maximum information. 

Accuracy of the calculation of Fourier coefficients. Our method of calculation of the Fourier 
coefficients is dependent on the knowledge of the rotation vector p, or an approximation p. For the rest 
of the chapter, meNis some fixed integer and M is some integer satisfying M > d + m(d + /3), where f3 is 
the Diophantine class of p. Let Ap be the error in the approximation of p, that is, 

Ap:=p-p. 

We are interested in knowing how the error in the calculation of Fourier coefficients depend on Ap. 
Let for every k e 7L d , fk(0) : = e l2nh0 . Then every periodic function g has the complex Fourier series 
representation 

g(6) = E a k f k , where a k = f g(0)f- k (0) 

k€ 1 d Jl d 

Therefore, a k can be approximated as 


a k = WB N p p (g(0)U(0)) = y 1 g(np)e- i2 ™ k ?w n>N . 

71=0 


Corollary 3.1 If N\\k ■ Ap|| << 1, then for each me N, there exists a constant C g , w , m > 0 such that 

\a k (g) - a k (g)\ 


MfiOl 


<nN\k-Ap\ + C g ^ m N~ m . 


Proof We begin by obtaining bounds on the error in calculating the Fourier coefficients of pure exponen¬ 
tials. Note that a k (fi ) equals 1 if l = k, equals 0 if l + k. We will estimate the error \a k (fi) -a k (f) \ in each 
case. 


MA) = WB^(/*(O/_ fc (0))| = E e~ i2 ™ k - A ?w mtN . 

71=0 

Then if N\\k • Ap\\ « 1, then using the approximation e ia - 1 » ia ) one can say that 

|^fc(//c) — ^k^fk) | — | ttk^fk) ~ 1| 

N -1 

= y( e -^.Ap_i 
77=0 

N -1 

/ \„-i27Tnk-Ap i \ * 

s L \ e - t\U!m,N 

77=0 

N-l 

< y 27 ink ■ \Ap\w m}N 

77=0 

= 7fN\k ■ Ap\. 
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Since p is Diophantine with Diophantine class ft > 0, for every a e M, there exists a constant C f) > 0 such 
that for all l e Z d with ||/|| sufficiently big, one has 

| exp(i27r(7 • p - a) - 1)| = | exp(f27r(7 ■ p - a) - exp(f27ro;))| > C'p||/|| _( ' <i+/3 ' > (24) 

By Eq. 14 from [3], for every l + k there exists a constant C Wtm > 0 depending on w and m such that for 
every «eR 

N-l 

Y^ exp(-i27rna)w m: N 

71=0 

Combining Eqs. [24] and [25j we get 

\a k (fi) ~ a k (fi)\ = \a k (fi)\ 

jv-i 

= exp(-f27rn(/ -p- k-p))w m , N 

71=0 

< C W}m N~ m \ exp(f27r(/ • p - k ■ p)) - l| _m 
<C w , m C p N- m \\l\\ m ^ 

Since / e C M , there exists C 9t M > 0 depending on g and M such that for each 11 0, |cq| < C'/ ! m|KI| _ jw - 
Therefore, in a manner similar to the derivation of Eq. 15 from [3], we can write 


< C w . m N~ m e 


-m I A2ixol 


-i|- 


(25) 


\a k (g ) - a k (g)\ = E Ma k (fi) - a k (f t )) 

le 7L d 

< \a k (a k (f k )-a k (f k ))\ + 


T\ai{a k (fi)-a k (fi)) 


< \a k \-nN\k • Ap| + E | ai(a k (fi) - a k (f t ))\ 

l±k 

< \aftirN\k ■ Ap\ + C w m CpN~ m Yi K|||Z|| m(d+/3) 

l±k 

< \a k \nN\k • Ap| + C w , m C p C E \\l\\- M \\l\\ m ^ 

l±k 

= \a k \nN\k ■ Ap\ + C w , m C p C fM N~ m E \\l\\-(M-m(d + p)) 


Since M > d + m(d + ft), the sum E \\l\\-( M - m ( d+ P)) < oo. The statement of the corollary now follows. 

l±k 


I 
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Figure 7: Two-dimensional torus map. Fig. (a) shows the first 10 4 iterates of an orbit for the two- 
dimensional quasiperiodic toms map. The orbit appears to be dense, indicating quasiperiodicity. This is 
confirmed by computing the two Lyapunov exponents of the orbit using Birkhoff averaging and finding 
them to be 0. The convergence of this computation for one of the two Lyapunov exponents is shown in 
blue in (b). The highest to lowest curves show the convergence rates resp. for the first three weighting 
functions given in (14). Fig. (c) shows the convergence rate for the first rotation number for the four 
different weighting functions. 
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Figure 8: Rate of convergence for different weight functions. Left is a plot of the three non-constant 
weighting functions from (14) (quad blue, sin 2 yellow, exp red). Since only the shape matters, they have 
been rescaled so that each has a peak of approximately one. For a given w and a given number of iterates 
N, the rotation number p approximation is calculated for Bi of the restricted three-body problem, the 
error of the calculation is the difference \p - p\. The figure shows the convergence rate of this error as 
a function of N. The exponential weight function is seen to be the best method. The error cannot be 
reduced below 10 -32 because that is the limit of quadruple precision. 




Figure 9: Conjugacy for the torus. Fig. (a) depicts the reconstruction of the periodic part g (see [ 5 ]) 
of the first component of the conjugacy function for the torus map. The second conjugacy function is 
similar but not depicted here. Fig. (b) shows the decay of the Fourier coefficients for this component of 
the conjugacy function on the log-linear scale. 
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